suppressPackageStartupMessages({
library(piecewiseSEM)
library(magrittr)
library(tidyverse)
library(stargazer)
library(broom)
})
storks <- read.csv("Data/storks.csv")
About the course
This course was part of the pre-conference workshops at the International Congress for Conservation Biology, in Cartagena, Colombia. The course was given by Dr. Achaz von Hardenberg @achazhardenberg.
Subjects on causal inference, path analalysis and bayesian inference were covered.
Storks deliver babies (p = 0.0079)
The scientist’s mantra: “Correlation does not imply causation”
Plot
ggplot(data = storks, mapping=aes(x = Storks, y = Birth)) +
geom_point() +
geom_smooth(method = "lm", color = "black") +
theme_bw() +
labs(x = "Number of breeding stork pairs", y = "Human birth rate (thousands/year)")

Linear model
lm(Birth ~ Storks, data = storks) %>%
stargazer(single.row = T, type = "html", dep.var.labels = "Human birth rate *thousands / year)")
|
|
|
|
Dependent variable:
|
|
|
|
|
|
Human birth rate *thousands / year)
|
|
|
|
Storks
|
0.029*** (0.009)
|
|
Constant
|
225.029** (93.561)
|
|
|
|
Observations
|
17
|
|
R2
|
0.385
|
|
Adjusted R2
|
0.344
|
|
Residual Std. Error
|
332.185 (df = 15)
|
|
F Statistic
|
9.380*** (df = 1; 15)
|
|
|
|
Note:
|
p<0.1; p<0.05; p<0.01
|
Other variables
Other variables appear to be correlated, so we can inspect their correlation and include relevant ones in the model.
library(corrplot)
corrplot(cor(storks[,-1]), method="ellipse", type="lower")

lm1 <- lm(Birth ~ Storks, storks)
lm2 <- lm(Area ~ Storks, storks)
lm3 <- lm(Birth ~ Area, storks)
lm4 <- lm(Birth ~ Area + Storks, storks)
stargazer(lm1, lm2, lm3, lm4, type ="html", single.row = T, add.lines = list(AIC = c("AIC", formatC(AIC(lm1, lm2, lm3, lm4)[,2], digits = 2, format = "f"))))
number of rows of result is not a multiple of vector length (arg 2)
|
|
|
|
Dependent variable:
|
|
|
|
|
|
Birth
|
Area
|
Birth
|
|
|
(1)
|
(2)
|
(3)
|
(4)
|
|
|
|
Storks
|
0.029*** (0.009)
|
14.400** (5.231)
|
|
0.006 (0.006)
|
|
Area
|
|
|
0.002*** (0.0002)
|
0.002*** (0.0002)
|
|
Constant
|
225.029** (93.561)
|
146,817.600** (52,057.720)
|
-7.775 (56.938)
|
-7.412 (56.702)
|
|
|
|
AIC
|
249.51
|
464.44
|
225.39
|
226.08
|
|
Observations
|
17
|
17
|
17
|
17
|
|
R2
|
0.385
|
0.336
|
0.851
|
0.862
|
|
Adjusted R2
|
0.344
|
0.291
|
0.841
|
0.842
|
|
Residual Std. Error
|
332.185 (df = 15)
|
184,830.100 (df = 15)
|
163.422 (df = 15)
|
162.744 (df = 14)
|
|
F Statistic
|
9.380*** (df = 1; 15)
|
7.578** (df = 1; 15)
|
85.731*** (df = 1; 15)
|
43.787*** (df = 2; 14)
|
|
|
|
Note:
|
p<0.1; p<0.05; p<0.01
|
After including Area, for example, we observe that the effect of Storks significantly decreases.
Correlation imples and unresolved causal structure
Once we have identified processes, we can start to imply causation:
Causality always imples a completely resolved correlation structure
Path analysis
Limitations of commonly employed statistical methods (i.e. multiple linear regression):
“Path analysis is an extension of multiple regression, and was developed to overcome these limitations”
In path analysis, we call each variable a vertex (or vertices, in plural) and their connections are called edges. Directed graphs include directions (i.e. the edges are arrow that define direction). In Indirect graphs, edges only connect vertices, but do not specify the direction of influence.
A full explanation on the grammar of path analysis is found in Achaz’s materials
From his model:
DiagrammeR::grViz("
digraph boxes_and_circles{
# Define nodes
node [shape = square
penwidth = 2,
fontsize = 24]
A; B; C; D; E; F
# Add edge statements
#Advaance classes
A -> B
B -> C
B -> D
E -> C
E -> D
F -> E
}
", height = 500)
A is:
B is:
C and D are:
E is:
F is:
In the case above, B is an active vertex, becauser it is bot an effect (of A) and a cause (of C or D). Similarly, D is an inactive vertez, which is only an effect of 1 or more other vertices, but not a cause of anything. Often, we try to understand relationships between D and A, leaving out B.
Exercise
Drawa causal graph with 6 variables:
DiagrammeR::grViz("
digraph boxes_and_circles{
# Define nodes
node [shape = square
penwidth = 2,
fontsize = 24]
A; B; C; D; E; F
# Add edge statements
#Advaance classes
A -> B
B -> C
C -> D
B -> D
D -> F
E -> F
}
", height = 500)
To start transforming this into statistical models, we perform d-separation:
- First, we identify all pairs of not adjacent vertices (i.e. not connected by a direct arrow)
d_separ <- data.frame(Pairs = c("C,D","A,C", "A,D", "A,E", "A,F", "B,E","B,F", "C,F","D,F"))
d_separ
- Then, we identify the parents of any of the vertices of each pair
d_separ <- data.frame(Pairs = c("C,D","A,C", "A,D", "A,E", "A,F", "B,E","B,F", "C,F","D,F"),
Parents = c("B,E", "B", "B,E", "F", "None", "A, F", "A", "B,E", "B,E"))
d_separ
- The D-separation statements are therefore:
d_separ %>%
mutate(D_separation = paste0("(", Pairs, ")", ", {", Parents, "}"))
The number of elements in the basis test should be given by
\[\frac{V!}{2(V-2)!}-A\]
In R code, we can define a function for this as:
numbers <- function(V, A){
n <- (factorial(V)/(2*factorial((V-2))))-A
return(n)
}
Where:
- V is the number of vertices
-A is the number of arrows
- Test probabilistic conditionals of every d-separation statement
d_separ %>%
mutate(D_separation = paste0("(", Pairs, ")", ", {", Parents, "}"),
Claims = c("D ~ B + E + C", "C ~ B + A", "D ~ B + E + A", "E ~ F + A", "F ~ A", "E ~ A + F + B", "F ~ B + A", "F ~ B + E + C", "F ~ B + E + D"))
We then use Fisher’s C test to test the composite probability of the whole set.
Re-visit storks and births
With a model of the sape:
DiagrammeR::grViz("
digraph boxes_and_circles{
# Define nodes
node [shape = square
penwidth = 2,
fontsize = 24]
A; S; B; I
# Add edge statements
#Advaance classes
A -> S
A -> B
B -> I
}
", height = 500)
Where:
A = Area
S = Storks
B = Birth rate
I = Inhabitants
Using the formula above, we know that we need to identify 3 pairs.
Using the steps above, we must identify all non-adjacent pairs
storke_pairs <- data.frame(Pairs = c("A, I", "S, B", "S,I"))
storke_pairs
Identify all parents of each pair
storke_pairs <- data.frame(Pairs = c("A, I", "S, B", "S,I"),
Parents = c("B", "A", "A, B"))
storke_pairs
The d-separation statements are then:
storke_pairs %>%
mutate(D_separation = paste0("(", Pairs, ")", ", {", Parents, "}"))
Test the probabilistic conditional dependence of every d-statement with the following claims:
storke_pairs %>%
mutate(D_separation = paste0("(", Pairs, ")", ", {", Parents, "}"),
Claims = c("I ~ B + A", "B ~ A + S", "I ~ A + B + I"))
Evaluate the models:
colnames(storks) <- c("Country", "A", "S", "I", "B")
path1 <- lm(I ~ B + A, storks)
path2 <- lm(B ~ A + S, storks)
path3 <- lm(I ~ A + B + S, storks)
stargazer(path1, path2, path3, type = "html", report = c("vcp"), single.row = T)
number of rows of result is not a multiple of vector length (arg 2)
|
|
|
|
Dependent variable:
|
|
|
|
|
|
I
|
B
|
I
|
|
|
(1)
|
(2)
|
(3)
|
|
|
|
B
|
0.039
|
|
0.049
|
|
|
p = 0.079
|
|
p = 0.032
|
|
A
|
0.00002
|
0.002
|
0.00002
|
|
|
p = 0.624
|
p = 0.00001
|
p = 0.575
|
|
S
|
|
0.006
|
-0.001
|
|
|
|
p = 0.307
|
p = 0.111
|
|
Constant
|
6.744
|
-7.412
|
6.771
|
|
|
p = 0.162
|
p = 0.898
|
p = 0.138
|
|
|
|
Observations
|
17
|
17
|
17
|
|
R2
|
0.729
|
0.862
|
0.779
|
|
Adjusted R2
|
0.691
|
0.842
|
0.728
|
|
Residual Std. Error
|
13.084 (df = 14)
|
162.744 (df = 14)
|
12.264 (df = 13)
|
|
F Statistic
|
18.872*** (df = 2; 14)
|
43.787*** (df = 2; 14)
|
15.297*** (df = 3; 13)
|
|
|
|
Note:
|
p<0.1; p<0.05; p<0.01
|
Using Fisher’s C:
\[C = -2\sum_{i = 1}^{k{}} log(p_i)\] Where:
- p is the rpobability of the term for wich we are testing independence (i.e. the other vertix from the non-adjancet verrix pairs)
C follows a \(\chi^2\) distribution, with 2k degrees of freedom: 1-pchisq(C,2*k)
extract_p <- function(model, coeff){
tidy(model) %>%
filter(term == coeff) %>%
select(p.value)
}
ppath1 <- extract_p(path1, "A")
ppath2 <- extract_p(path2, "S")
ppath3 <- extract_p(path3, "S")
C <- rbind(ppath1, ppath2, ppath3) %>%
mutate(logs = log(p.value)) %$%
sum(logs)
C <- -2*C
test <- 1-pchisq(C,2*3)
Thus, the probability is p = 0.2597204.
Some R-packages for all this
piecewiseSEM
Here, we only specify the dependent paths as within an lm() call. For example, in this model, I is caused by B, so we call lm(B~I). We include each lm call inside a list, and pass that to sem.fit, which will identify the missing paths, and do all the calculations.
list(
lm(S ~ A, storks),
lm(B ~ A, storks),
lm(I ~ B, storks)) %>%
sem.fit(data = storks, .progressBar = F, conditional = T) %>%
knitr::kable()
| I ~ A + B |
0e+00 |
0.0000 |
14 |
0.5024 |
0.6232 |
|
| B ~ S + A |
6e-03 |
0.0057 |
14 |
1.0609 |
0.3067 |
|
| I ~ S + A + B |
-8e-04 |
0.0004 |
13 |
-1.7128 |
0.1105 |
|
LS0tDQp0aXRsZTogIkNhdXNhbCBJbmZlcmVuY2UgZm9yIENvbnNlcnZhdGlvbiBCaW9sb2dpc3RzIg0KYXV0aG9yOiAiSnVhbiBDYXJsb3MgVmlsbGFzZcOxb3ItRGVyYmV6Ig0KZGF0ZTogIjIyIGRlIGp1bGlvIGRlIDIwMTciDQpvdXRwdXQ6IA0KICBodG1sX25vdGVib29rOg0KICAgIGNvZGVfZm9sZGluZzogaGlkZQ0KICAgIGZpZ19jYXB0aW9uOiB5ZXMNCiAgICB0b2M6IHllcw0KICAgIHRvY19jb2xsYXBzZTogbm8NCiAgICB0b2NfZmxvYXQ6IHllcw0KLS0tDQoNCmBgYHtyIHNldHVwLCBpbmNsdWRlPUZBTFNFfQ0Ka25pdHI6Om9wdHNfY2h1bmskc2V0KGVjaG8gPSBUUlVFKQ0KYGBgDQoNCmBgYHtyIGxvYWQgcGFja2FnZXN9DQpzdXBwcmVzc1BhY2thZ2VTdGFydHVwTWVzc2FnZXMoew0KICBsaWJyYXJ5KHBpZWNld2lzZVNFTSkNCiAgbGlicmFyeShtYWdyaXR0cikNCiAgbGlicmFyeSh0aWR5dmVyc2UpDQogIGxpYnJhcnkoc3RhcmdhemVyKQ0KICBsaWJyYXJ5KGJyb29tKQ0KfSkNCmBgYA0KDQpgYGB7cn0NCnN0b3JrcyA8LSByZWFkLmNzdigiRGF0YS9zdG9ya3MuY3N2IikNCmBgYA0KDQojIEFib3V0IHRoZSBjb3Vyc2UNCg0KVGhpcyBjb3Vyc2Ugd2FzIHBhcnQgb2YgdGhlIHByZS1jb25mZXJlbmNlIHdvcmtzaG9wcyBhdCB0aGUgSW50ZXJuYXRpb25hbCBDb25ncmVzcyBmb3IgQ29uc2VydmF0aW9uIEJpb2xvZ3ksIGluIENhcnRhZ2VuYSwgQ29sb21iaWEuIFRoZSBjb3Vyc2Ugd2FzIGdpdmVuIGJ5IERyLiBBY2hheiB2b24gSGFyZGVuYmVyZyBbXEBhY2hhemhhcmRlbmJlcmddKGh0dHBzOi8vdHdpdHRlci5jb20vYWNoYXpoYXJkZW5iZXJnKS4NCg0KU3ViamVjdHMgb24gY2F1c2FsIGluZmVyZW5jZSwgcGF0aCBhbmFsYWx5c2lzIGFuZCBiYXllc2lhbiBpbmZlcmVuY2Ugd2VyZSBjb3ZlcmVkLg0KDQojIFN0b3JrcyBkZWxpdmVyIGJhYmllcyAocCA9IDAuMDA3OSkNCg0KPiBUaGUgc2NpZW50aXN0J3MgbWFudHJhOiAqIkNvcnJlbGF0aW9uIGRvZXMgbm90IGltcGx5IGNhdXNhdGlvbiIqDQoNCiMjIFBsb3QNCg0KYGBge3J9DQpnZ3Bsb3QoZGF0YSA9IHN0b3JrcywgbWFwcGluZz1hZXMoeCA9IFN0b3JrcywgeSA9IEJpcnRoKSkgKw0KICBnZW9tX3BvaW50KCkgKw0KICBnZW9tX3Ntb290aChtZXRob2QgPSAibG0iLCBjb2xvciA9ICJibGFjayIpICsNCiAgdGhlbWVfYncoKSArDQogIGxhYnMoeCA9ICJOdW1iZXIgb2YgYnJlZWRpbmcgc3RvcmsgcGFpcnMiLCB5ID0gIkh1bWFuIGJpcnRoIHJhdGUgKHRob3VzYW5kcy95ZWFyKSIpDQpgYGANCg0KIyMgTGluZWFyIG1vZGVsDQoNCmBgYHtyLCByZXN1bHRzID0gImFzaXMifQ0KbG0oQmlydGggfiBTdG9ya3MsIGRhdGEgPSBzdG9ya3MpICU+JSANCiAgc3RhcmdhemVyKHNpbmdsZS5yb3cgPSBULCB0eXBlID0gImh0bWwiLCBkZXAudmFyLmxhYmVscyA9ICJIdW1hbiBiaXJ0aCByYXRlICp0aG91c2FuZHMgLyB5ZWFyKSIpDQpgYGANCg0KIyMgT3RoZXIgdmFyaWFibGVzDQoNCk90aGVyIHZhcmlhYmxlcyBhcHBlYXIgdG8gYmUgY29ycmVsYXRlZCwgc28gd2UgY2FuIGluc3BlY3QgdGhlaXIgY29ycmVsYXRpb24gYW5kIGluY2x1ZGUgcmVsZXZhbnQgb25lcyBpbiB0aGUgbW9kZWwuDQoNCmBgYHtyfQ0KbGlicmFyeShjb3JycGxvdCkNCg0KY29ycnBsb3QoY29yKHN0b3Jrc1ssLTFdKSwgbWV0aG9kPSJlbGxpcHNlIiwgdHlwZT0ibG93ZXIiKQ0KYGBgDQoNCmBgYHtyLCByZXN1bHRzID0gImFzaXMiLCBtZXNzYWdlID0gRn0NCmxtMSA8LSBsbShCaXJ0aCB+IFN0b3Jrcywgc3RvcmtzKQ0KbG0yIDwtIGxtKEFyZWEgfiBTdG9ya3MsIHN0b3JrcykNCmxtMyA8LSBsbShCaXJ0aCB+IEFyZWEsIHN0b3JrcykNCmxtNCA8LSBsbShCaXJ0aCB+IEFyZWEgKyBTdG9ya3MsIHN0b3JrcykNCg0Kc3RhcmdhemVyKGxtMSwgbG0yLCBsbTMsIGxtNCwgdHlwZSA9Imh0bWwiLCBzaW5nbGUucm93ID0gVCwgYWRkLmxpbmVzID0gbGlzdChBSUMgPSBjKCJBSUMiLCBmb3JtYXRDKEFJQyhsbTEsIGxtMiwgbG0zLCBsbTQpWywyXSwgZGlnaXRzID0gMiwgZm9ybWF0ID0gImYiKSkpKQ0KDQpgYGANCg0KQWZ0ZXIgaW5jbHVkaW5nIGBBcmVhYCwgZm9yIGV4YW1wbGUsIHdlIG9ic2VydmUgdGhhdCB0aGUgZWZmZWN0IG9mIGBTdG9ya3NgIHNpZ25pZmljYW50bHkgZGVjcmVhc2VzLg0KDQo+IENvcnJlbGF0aW9uIGltcGxlcyBhbmQgdW5yZXNvbHZlZCBjYXVzYWwgc3RydWN0dXJlDQoNCk9uY2Ugd2UgaGF2ZSBpZGVudGlmaWVkIHByb2Nlc3Nlcywgd2UgY2FuIHN0YXJ0IHRvIGltcGx5IGNhdXNhdGlvbjoNCg0KPiBDYXVzYWxpdHkgYWx3YXlzIGltcGxlcyBhIGNvbXBsZXRlbHkgcmVzb2x2ZWQgY29ycmVsYXRpb24gc3RydWN0dXJlDQoNCiMgUGF0aCBhbmFseXNpcw0KDQpMaW1pdGF0aW9ucyBvZiBjb21tb25seSBlbXBsb3llZCBzdGF0aXN0aWNhbCBtZXRob2RzICgqaS5lLiogbXVsdGlwbGUgbGluZWFyIHJlZ3Jlc3Npb24pOg0KDQotIENhbiBvbmx5IGFuYWx5emUgb25lIGRlcGVuZGVudCB2YXJpYWJsZSBhdCBhIHRpbWUNCg0KLSBBIHBhcnRpY3VsYXIgdmFyaWFibGUgY2FuIGVpdGhlciBiZSBhIHByZWRpY3RvciBvciBhIHJlc3BvbnNlDQoNCioiUGF0aCBhbmFseXNpcyBpcyBhbiBleHRlbnNpb24gb2YgbXVsdGlwbGUgcmVncmVzc2lvbiwgYW5kIHdhcyBkZXZlbG9wZWQgdG8gb3ZlcmNvbWUgdGhlc2UgbGltaXRhdGlvbnMiKg0KDQpJbiBwYXRoIGFuYWx5c2lzLCB3ZSBjYWxsIGVhY2ggdmFyaWFibGUgYSB2ZXJ0ZXggKG9yIHZlcnRpY2VzLCBpbiBwbHVyYWwpIGFuZCB0aGVpciBjb25uZWN0aW9ucyBhcmUgY2FsbGVkIGVkZ2VzLiBEaXJlY3RlZCBncmFwaHMgaW5jbHVkZSBkaXJlY3Rpb25zICgqaS5lLiogdGhlIGVkZ2VzIGFyZSBhcnJvdyB0aGF0IGRlZmluZSBkaXJlY3Rpb24pLiBJbiBJbmRpcmVjdCBncmFwaHMsIGVkZ2VzIG9ubHkgY29ubmVjdCB2ZXJ0aWNlcywgYnV0IGRvIG5vdCBzcGVjaWZ5IHRoZSBkaXJlY3Rpb24gb2YgaW5mbHVlbmNlLg0KDQpBIGZ1bGwgZXhwbGFuYXRpb24gb24gdGhlIGdyYW1tYXIgb2YgcGF0aCBhbmFseXNpcyBpcyBmb3VuZCBpbiBbQWNoYXoncyBtYXRlcmlhbHNdKCJDYXVzYWwgaW5mZXJlbmNlIENhcnRhZ2VuYSAucGRmIikNCg0KRnJvbSBoaXMgbW9kZWw6DQoNCmBgYHtyfQ0KRGlhZ3JhbW1lUjo6Z3JWaXooIg0KICAgICAgZGlncmFwaCBib3hlc19hbmRfY2lyY2xlc3sNCg0KIyBEZWZpbmUgbm9kZXMNCm5vZGUgW3NoYXBlID0gc3F1YXJlDQogICAgICBwZW53aWR0aCA9IDIsDQogICAgICBmb250c2l6ZSA9IDI0XQ0KDQpBOyBCOyBDOyBEOyBFOyBGDQoNCiMgQWRkIGVkZ2Ugc3RhdGVtZW50cw0KI0FkdmFhbmNlIGNsYXNzZXMNCkEgLT4gQg0KQiAtPiBDDQpCIC0+IEQNCkUgLT4gQw0KRSAtPiBEDQpGIC0+IEUNCn0NCiAgICAgICIsIGhlaWdodCA9IDUwMCkNCmBgYA0KDQpBIGlzOg0KDQotIGRpcmVjdCBjYXVzZSBvZiBCDQoNCi0gaW5kaXJlY3QgY2F1c2Ugb2YgQyBhbmQgRg0KDQpCIGlzOg0KDQotIGNhdXNlZCBieSBBDQoNCkMgYW5kIEQgYXJlOg0KDQotIGluZGlyZWN0bHkgaW5kZXBlbmRlbnQgb24gQSBhbmQgRg0KDQotIGRpcmVjdGx5IGRlcGVuZGVudCBvbiBCIGFuZCBFDQoNCkUgaXM6DQoNCi0gZGlyZWN0bHkgZGVwZW5kZW50IG9uIEYNCg0KRiBpczoNCg0KLSBkaXJlY3QgY2F1c2Ugb2YgRQ0KDQotIEluZGlyZWN0IGNhdXNlIG9mIEMgYW5kIEYNCg0KSW4gdGhlIGNhc2UgYWJvdmUsIEIgaXMgYW4gYWN0aXZlIHZlcnRleCwgYmVjYXVzZXIgaXQgaXMgYm90IGFuIGVmZmVjdCAob2YgQSkgYW5kIGEgY2F1c2UgKG9mIEMgb3IgRCkuIFNpbWlsYXJseSwgRCBpcyBhbiBpbmFjdGl2ZSB2ZXJ0ZXosIHdoaWNoIGlzIG9ubHkgYW4gZWZmZWN0IG9mIDEgb3IgbW9yZSBvdGhlciB2ZXJ0aWNlcywgYnV0IG5vdCBhIGNhdXNlIG9mIGFueXRoaW5nLiBPZnRlbiwgd2UgdHJ5IHRvIHVuZGVyc3RhbmQgcmVsYXRpb25zaGlwcyBiZXR3ZWVuIEQgYW5kIEEsIGxlYXZpbmcgb3V0IEIuDQoNCiMjIyBFeGVyY2lzZQ0KDQpEcmF3YSBjYXVzYWwgZ3JhcGggd2l0aCA2IHZhcmlhYmxlczoNCg0KLSAoZnJvbUF0b0Ypd2hlcmU6DQoNCi0gQSBpcyBkaXJlY3QgY2F1c2Ugb2YgQg0KDQotIEMgaXMgYSBjYXVzYWwgY2hpbGQgb2YgQg0KDQotIEQgaXMgZGlyZWN0bHkgZGVwZW5kZW50IG9uIEIgYW5kIEMNCg0KLSBGIGlzIGRpcmVjdGx5IGRlcGVuZGVudCBvbiBEIGFuZCBFDQoNCmBgYHtyfQ0KRGlhZ3JhbW1lUjo6Z3JWaXooIg0KICAgICAgZGlncmFwaCBib3hlc19hbmRfY2lyY2xlc3sNCg0KIyBEZWZpbmUgbm9kZXMNCm5vZGUgW3NoYXBlID0gc3F1YXJlDQogICAgICBwZW53aWR0aCA9IDIsDQogICAgICBmb250c2l6ZSA9IDI0XQ0KDQpBOyBCOyBDOyBEOyBFOyBGDQoNCiMgQWRkIGVkZ2Ugc3RhdGVtZW50cw0KI0FkdmFhbmNlIGNsYXNzZXMNCkEgLT4gQg0KQiAtPiBDDQpDIC0+IEQNCkIgLT4gRA0KRCAtPiBGDQpFIC0+IEYNCg0KfQ0KICAgICAgIiwgaGVpZ2h0ID0gNTAwKQ0KYGBgDQoNClRvIHN0YXJ0IHRyYW5zZm9ybWluZyB0aGlzIGludG8gc3RhdGlzdGljYWwgbW9kZWxzLCB3ZSBwZXJmb3JtIGQtc2VwYXJhdGlvbjoNCg0KMS4gRmlyc3QsIHdlIGlkZW50aWZ5IGFsbCBwYWlycyBvZiBub3QgYWRqYWNlbnQgdmVydGljZXMgKCppLmUuIG5vdCBjb25uZWN0ZWQgYnkgYSBkaXJlY3QgYXJyb3cqKQ0KDQpgYGB7cn0NCmRfc2VwYXIgPC0gZGF0YS5mcmFtZShQYWlycyA9IGMoIkMsRCIsIkEsQyIsICJBLEQiLCAiQSxFIiwgIkEsRiIsICJCLEUiLCJCLEYiLCAiQyxGIiwiRCxGIikpDQoNCmRfc2VwYXINCmBgYA0KDQoyLiBUaGVuLCB3ZSBpZGVudGlmeSB0aGUgcGFyZW50cyBvZiBhbnkgb2YgdGhlIHZlcnRpY2VzIG9mIGVhY2ggcGFpcg0KDQpgYGB7cn0NCg0KZF9zZXBhciA8LSBkYXRhLmZyYW1lKFBhaXJzID0gYygiQyxEIiwiQSxDIiwgIkEsRCIsICJBLEUiLCAiQSxGIiwgIkIsRSIsIkIsRiIsICJDLEYiLCJELEYiKSwNCiAgICAgICAgICAgICAgICAgICAgICBQYXJlbnRzID0gYygiQixFIiwgIkIiLCAiQixFIiwgIkYiLCAiTm9uZSIsICJBLCBGIiwgIkEiLCAiQixFIiwgIkIsRSIpKQ0KDQpkX3NlcGFyDQpgYGANCg0KMy4gVGhlIEQtc2VwYXJhdGlvbiBzdGF0ZW1lbnRzIGFyZSB0aGVyZWZvcmU6DQoNCmBgYHtyfQ0KZF9zZXBhciAlPiUgDQogIG11dGF0ZShEX3NlcGFyYXRpb24gPSBwYXN0ZTAoIigiLCBQYWlycywgIikiLCAiLCB7IiwgUGFyZW50cywgIn0iKSkNCmBgYA0KDQpUaGUgbnVtYmVyIG9mIGVsZW1lbnRzIGluIHRoZSBiYXNpcyB0ZXN0IHNob3VsZCBiZSBnaXZlbiBieQ0KDQokJFxmcmFje1YhfXsyKFYtMikhfS1BJCQNCg0KSW4gUiBjb2RlLCB3ZSBjYW4gZGVmaW5lIGEgZnVuY3Rpb24gZm9yIHRoaXMgYXM6DQoNCmBgYHtyfQ0KbnVtYmVycyA8LSBmdW5jdGlvbihWLCBBKXsNCiAgbiA8LSAoZmFjdG9yaWFsKFYpLygyKmZhY3RvcmlhbCgoVi0yKSkpKS1BDQogIHJldHVybihuKQ0KfQ0KYGBgDQoNCldoZXJlOg0KDQotIFYgaXMgdGhlIG51bWJlciBvZiB2ZXJ0aWNlcw0KDQotQSBpcyB0aGUgbnVtYmVyIG9mIGFycm93cw0KDQo0LiBUZXN0IHByb2JhYmlsaXN0aWMgY29uZGl0aW9uYWxzIG9mIGV2ZXJ5IGQtc2VwYXJhdGlvbiBzdGF0ZW1lbnQNCg0KYGBge3J9DQpkX3NlcGFyICU+JSANCiAgbXV0YXRlKERfc2VwYXJhdGlvbiA9IHBhc3RlMCgiKCIsIFBhaXJzLCAiKSIsICIsIHsiLCBQYXJlbnRzLCAifSIpLA0KICAgICAgICAgQ2xhaW1zID0gYygiRCB+IEIgKyBFICsgQyIsICJDIH4gQiArIEEiLCAiRCB+IEIgKyBFICsgQSIsICJFIH4gRiArIEEiLCAiRiB+IEEiLCAiRSB+IEEgKyBGICsgQiIsICJGIH4gQiArIEEiLCAiRiB+IEIgKyBFICsgQyIsICJGIH4gQiArIEUgKyBEIikpDQoNCmBgYA0KDQpXZSB0aGVuIHVzZSBGaXNoZXIncyBDIHRlc3QgdG8gdGVzdCB0aGUgY29tcG9zaXRlIHByb2JhYmlsaXR5IG9mIHRoZSB3aG9sZSBzZXQuDQogDQojIFJlLXZpc2l0IHN0b3JrcyBhbmQgYmlydGhzDQogDQpXaXRoIGEgbW9kZWwgb2YgdGhlIHNhcGU6DQogDQogDQogDQpgYGB7cn0NCkRpYWdyYW1tZVI6OmdyVml6KCINCiAgICAgIGRpZ3JhcGggYm94ZXNfYW5kX2NpcmNsZXN7DQoNCiMgRGVmaW5lIG5vZGVzDQpub2RlIFtzaGFwZSA9IHNxdWFyZQ0KICAgICAgcGVud2lkdGggPSAyLA0KICAgICAgZm9udHNpemUgPSAyNF0NCg0KQTsgUzsgQjsgSQ0KDQojIEFkZCBlZGdlIHN0YXRlbWVudHMNCiNBZHZhYW5jZSBjbGFzc2VzDQpBIC0+IFMNCkEgLT4gQg0KQiAtPiBJDQoNCn0NCiAgICAgICIsIGhlaWdodCA9IDUwMCkNCmBgYA0KDQpXaGVyZToNCg0KLSAgQSA9IEFyZWENCg0KLSBTID0gU3RvcmtzDQoNCi0gQiA9IEJpcnRoIHJhdGUNCg0KLSBJID0gSW5oYWJpdGFudHMNCg0KVXNpbmcgdGhlIGZvcm11bGEgYWJvdmUsIHdlIGtub3cgdGhhdCB3ZSBuZWVkIHRvIGlkZW50aWZ5IGByIG51bWJlcnMoNCwgMylgIHBhaXJzLg0KDQpVc2luZyB0aGUgc3RlcHMgYWJvdmUsIHdlIG11c3QgaWRlbnRpZnkgYWxsIG5vbi1hZGphY2VudCBwYWlycw0KDQpgYGB7cn0NCnN0b3JrZV9wYWlycyA8LSBkYXRhLmZyYW1lKFBhaXJzID0gYygiQSwgSSIsICJTLCBCIiwgIlMsSSIpKQ0KDQpzdG9ya2VfcGFpcnMNCmBgYA0KDQpJZGVudGlmeSBhbGwgcGFyZW50cyBvZiBlYWNoIHBhaXINCg0KYGBge3J9DQpzdG9ya2VfcGFpcnMgPC0gZGF0YS5mcmFtZShQYWlycyA9IGMoIkEsIEkiLCAiUywgQiIsICJTLEkiKSwNCiAgICAgICAgICAgICAgICAgICAgICAgICAgIFBhcmVudHMgPSBjKCJCIiwgIkEiLCAiQSwgQiIpKQ0KDQpzdG9ya2VfcGFpcnMNCmBgYA0KDQpUaGUgZC1zZXBhcmF0aW9uIHN0YXRlbWVudHMgYXJlIHRoZW46DQpgYGB7cn0NCnN0b3JrZV9wYWlycyAlPiUgDQogIG11dGF0ZShEX3NlcGFyYXRpb24gPSBwYXN0ZTAoIigiLCBQYWlycywgIikiLCAiLCB7IiwgUGFyZW50cywgIn0iKSkNCmBgYA0KDQpUZXN0IHRoZSBwcm9iYWJpbGlzdGljIGNvbmRpdGlvbmFsIGRlcGVuZGVuY2Ugb2YgZXZlcnkgZC1zdGF0ZW1lbnQgd2l0aCB0aGUgZm9sbG93aW5nIGNsYWltczoNCg0KYGBge3J9DQpzdG9ya2VfcGFpcnMgJT4lIA0KICBtdXRhdGUoRF9zZXBhcmF0aW9uID0gcGFzdGUwKCIoIiwgUGFpcnMsICIpIiwgIiwgeyIsIFBhcmVudHMsICJ9IiksDQogICAgICAgICBDbGFpbXMgPSBjKCJJIH4gQiArIEEiLCAiQiB+IEEgKyBTIiwgIkkgfiBBICsgQiArIEkiKSkNCmBgYA0KDQpFdmFsdWF0ZSB0aGUgbW9kZWxzOg0KDQpgYGB7ciwgcmVzdWx0cyA9ICJhc2lzIn0NCg0KY29sbmFtZXMoc3RvcmtzKSA8LSBjKCJDb3VudHJ5IiwgIkEiLCAiUyIsICJJIiwgIkIiKQ0KcGF0aDEgPC0gbG0oSSB+IEIgKyBBLCBzdG9ya3MpDQpwYXRoMiA8LSBsbShCIH4gQSArIFMsIHN0b3JrcykNCnBhdGgzIDwtIGxtKEkgfiBBICsgQiArIFMsIHN0b3JrcykNCg0Kc3RhcmdhemVyKHBhdGgxLCBwYXRoMiwgcGF0aDMsIHR5cGUgPSAiaHRtbCIsIHJlcG9ydCA9IGMoInZjcCIpLCBzaW5nbGUucm93ID0gVCkNCmBgYA0KDQpVc2luZyBGaXNoZXIncyBDOg0KDQokJEMgPSAtMlxzdW1fe2kgPSAxfV57a3t9fSBsb2cocF9pKSQkDQpXaGVyZToNCg0KLSAgcCBpcyB0aGUgcnBvYmFiaWxpdHkgb2YgdGhlIHRlcm0gZm9yIHdpY2ggd2UgYXJlIHRlc3RpbmcgaW5kZXBlbmRlbmNlICgqaS5lLiogdGhlIG90aGVyIHZlcnRpeCBmcm9tIHRoZSBub24tYWRqYW5jZXQgdmVycml4IHBhaXJzKQ0KDQpDIGZvbGxvd3MgYSAkXGNoaV4yJCBkaXN0cmlidXRpb24sIHdpdGggMmsgZGVncmVlcyBvZiBmcmVlZG9tOiBgMS1wY2hpc3EoQywyKmspYA0KDQpgYGB7cn0NCmV4dHJhY3RfcCA8LSBmdW5jdGlvbihtb2RlbCwgY29lZmYpew0KICB0aWR5KG1vZGVsKSAlPiUNCiAgICBmaWx0ZXIodGVybSA9PSBjb2VmZikgJT4lDQogICAgc2VsZWN0KHAudmFsdWUpDQp9DQoNCnBwYXRoMSA8LSBleHRyYWN0X3AocGF0aDEsICJBIikNCnBwYXRoMiA8LSBleHRyYWN0X3AocGF0aDIsICJTIikNCnBwYXRoMyA8LSBleHRyYWN0X3AocGF0aDMsICJTIikNCg0KQyA8LSByYmluZChwcGF0aDEsIHBwYXRoMiwgcHBhdGgzKSAlPiUgDQogIG11dGF0ZShsb2dzID0gbG9nKHAudmFsdWUpKSAlJCUNCiAgc3VtKGxvZ3MpDQoNCkMgPC0gIC0yKkMNCg0KdGVzdCA8LSAxLXBjaGlzcShDLDIqMykNCmBgYA0KDQoNClRodXMsIHRoZSBwcm9iYWJpbGl0eSBpcyAqcCA9ICogYHIgdGVzdGAuDQoNCiMgU29tZSBSLXBhY2thZ2VzIGZvciBhbGwgdGhpcw0KDQojIyBgcGllY2V3aXNlU0VNYA0KDQpIZXJlLCB3ZSBvbmx5IHNwZWNpZnkgdGhlIGRlcGVuZGVudCBwYXRocyBhcyB3aXRoaW4gYW4gYGxtKClgIGNhbGwuIEZvciBleGFtcGxlLCBpbiB0aGlzIG1vZGVsLCBJIGlzIGNhdXNlZCBieSBCLCBzbyB3ZSBjYWxsIGBsbShCfkkpYC4gV2UgaW5jbHVkZSBlYWNoIGxtIGNhbGwgaW5zaWRlIGEgbGlzdCwgYW5kIHBhc3MgdGhhdCB0byBgc2VtLmZpdGAsIHdoaWNoIHdpbGwgaWRlbnRpZnkgdGhlIG1pc3NpbmcgcGF0aHMsIGFuZCBkbyBhbGwgdGhlIGNhbGN1bGF0aW9ucy4NCg0KYGBge3J9DQpsaXN0KA0KICBsbShTIH4gQSwgc3RvcmtzKSwNCiAgbG0oQiB+IEEsIHN0b3JrcyksDQogIGxtKEkgfiBCLCBzdG9ya3MpKSAlPiUgDQogIHNlbS5maXQoZGF0YSA9IHN0b3JrcywgLnByb2dyZXNzQmFyID0gRiwgY29uZGl0aW9uYWwgPSBUKSAlPiUgDQogIGtuaXRyOjprYWJsZSgpDQpgYGANCg0K